Simulation of many-qubit quantum computation with matrix product states 
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Matrix product states provide a natural entanglement basis to represent a quantum register and 
operate quantum gates on it. This scheme can be materialized to simulate a quantum adiabatic 
algorithm solving hard instances of a NP-Complete problem. Errors inherent to truncations of the 
exact action of interacting gates are controlled by the size of the matrices in the representation. 
The property of finding the right solution for an instance and the expected value of the energy 
(cost function) are found to be remarkably robust against these errors. As a symbolic example, 
we simulate the algorithm solving a 100-qubit hard instance, that is, finding the correct product 
state out of ~ 10 30 possibilities. Accumulated statistics for up to 60 qubits seem to point at a 
sub-exponential growth of the average minimum time to solve hard instances with highly-truncated 
simulations of adiabatic quantum evolution. 

PACS numbers: 03.67.-a, 03.65.Ud, 03.67.Hk 
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A detailed understanding of a many-spin quantum sys- 
tem often requires its simulation on a classical computer. 
Such a possibility is limited to a small number of spins 
due to the exponential growth of the size of the Hilbcrt 
space. This is at the heart of the motivation to build a 
quantum computer 0. Using Standard present technol- 
ogy, a faithful simulation of a general Hamiltonian can 
be achieved for systems up to the order of 24 spins. 

Recent developments in representing quantum states 
and operating unitary evolution on them have refined the 
above common lore. The idea has evolved from accumu- 
lated knowledge on matrix product states (MPS, related 
to the density matrix renormalization group techniquc) 
2] and new insights from quantum information. Let us 
recali that a quantum state for an n-qubit system can be 
represented by the matrix product construction 
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where the indices i\, . . . , i n for each qubit range from 
to 1 (the qubits are placed in a chain) and oti, . . . , a n -\ 
are referred to as ancillae indices that range from 1 to a 
parameter we shall call Each matrix A^llaa at site a 
can be viewed as a projector from a pair of unphysical an- 
cillae to the physical degree of freedom that we associate 
to the computational basis. The success of MPS con- 
sists in changing the representation of the quantum state 
from the computational basis to a non-local one, closely 
attached to entanglement. To make this comment con- 
crete, let us note that the matrix representation of a state 
can be recovered via a chain of Schmidt dccompositions 
that separate a local system at a time, as made explicit 



by Vidal |3j. More specifically, A K a 
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Aa a being the Schmidt coefficients of the cut of the sys- 
tem between the a and a + 1 sides, and r( a ' being tensors 
for qubit a. The larger the entanglement is for different 



partitions of the system, the larger is the needed ancil- 
lae space, which corresponds to a higher rank \- MPS 
can handle simulations of various dynamics of spin chains 
with up to hundreds of spins because their little amount 
of entanglement can be represented with \ = 0(poly(n)) 
0- Q!| ■ A number of new developments have popped up 
from the bàsic MPS in the context of quantum informa- 
tion. In rcf . 3] , an efScient implementation of Hamilto- 
nian evolution was constructed for slightly entangled sys- 
tems. An explicit renormalization group transformation 
on quantum states was made explicit using MPS 5]. The 
rigid linear structure of MPS is being now abandoned in 
favor of the more general projected entangled-pair states 
(PEPS) that have been successfully applied to higher di- 
mensional systems 0. 

The natural question arises of whether MPS can be 
applied to simulate a quantum computer. The content of 
this paper is aimed to show that this is indeed possible 
and that we can handle large simulations with controlled 
aceuracy. As we shall describe, each time an entangling 
gate is operated on two neighboring qubits, the range of 
the connected ancillae index is doubled. This is the way 
interacting gates entangle the system. To keep the simu- 
lation under control, a (non-unique) truncation scheme is 
needed that stops the exponential growth of ancillae di- 
mensions. We expect this approximation scheme to fail 
whenever the inherently needed x is 0(2"). Nevcrthc- 
less, in some of these cases keeping \ = 0(poly(n)) in 
the simulation already gives reasonable approximations 
to the exact calculation, as we shall see. 

Our presentation will be made concrete by showing 
an MPS simulation of quantum computation in the case 
of adiabatic evolution for the NP-Complete Exact Cover 
satisfiability problem 0, • An instance of Exact Cover 
is defined by a set of m 3-bit clauses with satisfying as- 
signments 001, 010 or 100. The problem is defined as 
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deciding whether a given instance accepts a global satis- 
fying assignment of n bits. This satisfiability problem is 
NP-Complete. Classically hard instances of Exact Cover 
seem to appear at the so-called easy-hard-easy transition 
around m ~ .8ra ,9j. We have constructed such hard in- 
stances, with the additional property of having a uniquc 
satisfying assignment. The generation of hard instances 
is in itself a difncult problem for which we have developed 
specific algorithms, all of them based on the iterative ad- 
dition of random clauses that strictly decrease the num- 
ber of solutions of the instance until a single satisfying 
assignment is reached. 

The quantum algorithm for a given Exact Cover in- 
stance follows the adiabatic evolution of the ground state 
of a Hamiltonian (cost operator) defined by H{s) = (1 — 
s)Hq + sí/p, where the adiabatic parameter is s = t/T 
and t runs up to a total predetermined time T . We take 
the initial Hamiltonian to be Hq = X^ILi if (1 — a f ) where 
di stands for the number of clauses qubit i enters. The 
non-local problem Hamiltonian corresponds to the sum of 
clauses defined as Hp = X) c (i j k) ( z i + z i + z k ~ l) 2 where 
Zi = (1 — of)/2 has eigenvalues and 1, and c(i,j,k) 
stands for a clause involving qubits i, j and k. Exact 
simulations of quantum algorithms by adiabatic evolu- 
tion solving hard instances of satisfiability problems have 
been carried so far up to 30 qubits 0. The explosion 
of entanglement between random cuts in the quantum 
register was first analyzed in ref. The adiabatic 

evolution drives the system near a quantum phase tran- 
sition at s ~ .69 following universal sealing laws. Entropy 
for half-cuts of the register approximates on average the 
sealing law S ~ An, which almost saturates the maxi- 
mum S = n/2. This implies that the quantum algorithm 
cannot be simulated efficiently in a classical computer 
Ü|. Yet, the fact that entropy does not reach its al- 
lowed maximum suggests that an adequate handling of 
entanglement may provide a way to extend simulations 
far from naive limitations. 

Let us now turn to discuss the detailed way MPS can 
handle the simulation of the adiabatic evolution of Ex- 
act Cover. The simulation needs to follow a time evo- 
lution controlled by the s-dependent Hamiltonian. This 
continuous unitary time evolution can be discretized as 
follows: Ut,o = Ut,t-a ■ ■ ■ C^2A,aC^a,o where the incre- 
ment A = defines the discretization, M being a pos- 
itive integer. Our simulations indicate that we can take 
A = 0.125 while keeping sufRcient aceuracy (as com- 
pared to smaller A) in all of them. We have explicitely 
checked that simulations performed with A < 0.125 lead 
to equally-good discretizations of the continuous-time 
adiabatic algorithm, in the sense that the obtained re- 
sults do not practically differ from the ones calculated 
for A = 0.125. After l steps s = ± = ^ = ±, be- 
ing l — 0, . . . M. At any point in the evolution Trotter's 
formula to second order is used to divide the unitary op- 
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into elementary gates: í7(/+i)a,ía 



e iAH(s) ^ f e i§(i-s)H 0e iSsHp gif (í-s)-fío^ T _ We have ver- 

ified that we can maintain a faithful simulation with 
5 = A. The split of exponentials in Trotter's expansion 
is chosen so that Hq is separated from Hp. This brings 
the advantage that each piece of the evolution operator 
can be decomposed in mutually commuting elementary 
gates: 
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The adiabatic evolution is thus finally reduced to a series 
of one and two-qubit gates. The detailed way these gates 
operate on the MPS follows the original idea of ref. Q: 

1. A one-qubit gate acting on qubit a only involves an 
updating of that goes as follows: 

U^A^\i a )=A^üQy a ), (4) 
which corresponds to the local updating rule 

A ad - U i a i' a A aa ■ W 

This gate does not affect ancillae indices. Entanglement 
is unaffected as we are just performing local operations. 

As an example, consider the one-qubit gate = 
ei '* , cr^' being the usual Pauli matrix 
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acting on qubit a. Then, we have the following simple 
updating rule for A^: 
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2. A two-qubit gate involving contiguous qubits a and 
a + 1 follows a similar strategy. Let us define 
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Unlikc one-qubit gates, interacting gates do not preserve 
the product form of the tensors A. To reestablish the 
MPS structure we need to rewrite <d using a Schmidt 
decomposition. The procedure to follow is to computo 
the reduced density matrix from the cut of the system 
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between the a and a + 1 sites, which for the right side 
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, where we have made use 
(o-i) 



reads p 1 ^ 

of the x known Schmidt coefficients A^" - ^ for the cut 
between the a — 1 and the a sites. After the diagonaliza- 
tion of p using (ia) and (j'7) as composed indices, we di- 
rectly read from the eigenvalues the updated 2\ Schmidt 
coefncients A'l for this cut, and the updated matrices 

from the coefncients of the eigenvectors. Fi- 
nally, the new tensors for qubit a are easily calculated as 

A/(a+l)Í a +l G.Ía.Ía+1 



4/(0) 
^ a/3 
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Let us clarify this procedure with a simple example 
consider the quantum state of two qubits 



m = 100) 



(9) 



It is easy to verify that the above state is described by 
the following vàlues of the matrices A'· a > : 
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Notice that since the state is separable x = 1. 
point, let us apply the two-qubit gate 
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(10) 
At this 
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to the quantum state | 
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(|00> + |H>) 
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Since the resultant state is a maximaly entangled state 
of 2 qubits, we expect x to be bigger than 1. In order to 
evaluate the updated matrices A'^ a ' for qubits 1 and 2 
we compute the quantity defined in equation ||SJ|, which 
in our case turns out to be 



üü = e n = 



e 



01 



1 

V2 
1 

V2 



(13) 



The density matrix for qubit 2 (which in this case is 
equivalent to the density matrix for qubit 1) then reads 
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Since the above density matrix is already diagonal, it is 
clear that the updated Schmidt coefncients will be 



A'í 1} = A' 
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and, as expected, we see that x = 2 since entanglement 
has been created by the two-qubit gate. From the above 



expressions it is simple to get the value of the updated 
matrices A'^ for qubits 1 and 2: 
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3. Operations involving non-contiguous qubits (as in 
Exact Cover clauses) can be reduced to the case 2 us- 
ing SWAP operations, producing an overhead of 0(n) 
operations per clause. 

The exact simulation of a quantum computer is then 
completely defined. The running time of this algorithm 
scales as ~ Tnmx 3 - Efficiency depends on the way the 
growth of the ancillae space is handled. To keep the sim- 
ulation under control we define a truncation scheme of 
the exact simulation. We choose to use a local proce- 
dure, namely, we keep the first x terms out of the 2x in 
the Schmidt decomposition defined in the point 2 above. 
Only the terms that carry most of the entanglement in 
the decomposition are kept |3j. This reasonable trunca- 
tion carries an inherent -but always under control- loss 
of unitarity, since the sum of the retained squared eigen- 
values will not reach 1. As we shall see, larger x's allow 
for more faithful simulations. Alternatively, it would be 
possible to recast the whole enlarged state into its origi- 
nal size in an optimal way Q . While this second method 
is manifestly more precise, it carries an operational time 
overhead. It is then worth analyzing both techniques. In 
this paper we shall focus on the first one and leave the 
second for a separate publication. 

We have implemented a number of optimizations upon 
the above bàsic scheme. For any non-local gate there is 
an overhead of SWAP operations that damage the preci- 
sion of the computation. To minimize this effect, every 
three-qubit clause is operated as follows: we bring to- 
gether the three qubits with SWAPs of the left and right 
qubits keeping the central one fixed and, then, we oper- 
ate the two-qubit gates. Before returning the qubits to 
their original position we check if any of them is needed 
in the next gate. If so, we save whatever SWAP may be 
compensated between the two gates. Ordering of gates 
is also used to produce a saving of ~ 2/3 of the naive 
SWAPs. Diagonalization of the density matrix in the 
minimum allowed Hilbert space is used as well. A furthcr 
improvement is to keep a dynamical and local \, so that 
ancillae indices at the different partitions are allowed to 
take independent vàlues and grow up to site-dependent 
limits. This procedure, though, has shown essentially no 
improvement upon a fixed-x strategy. 

Let us now focus on the results. We first simulate 
the adiabatic algorithm with the requirement that the 
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Figure 1: (Color online) Computation of the absolute error 
(compared to the \ = 40 case) of the expected value of the 
energy (cost function) along the adiabatic evolution for a typ- 
ical instance with 30 qubits and 24 clauses for T = 100 as \ 
increases. Note the increasing precision with larger \ as s 
approaches the phase transition from the left-hand-side. In 
the inset, the absolute cost function is plotted. A similar be- 
havior is also obtained for other instances, all of them getting 
the exact solution at the end of the computation. 




Figure 2: (Color online) Norm in the register as a function 
of x m logarithmic scale, for instances of 14, 18, 22 and 30 
qubits. 



right solution is found for a typical instance of n = 30 
qubits with m = 24 clauses and T = 100. Along the 
evolution we compute the expected value of the energy 
(cost function) of the system, which can be calculated in 
0(n poly(x)) time. This is shown in Fig. ^ The system 
remains remarkably close to the instantaneous ground- 
state cost function along the approximated evolution. 
The error in the cost function is minimizcd as x increases. 
It is noteworthy to observe how the error in the simula- 
tion of the adiabatic algorithm peakes at the phase tran- 
sition point. We have also checked that it is precisely 
at this point where each qubit makes a decision towards 



Figure 3: (Color online) An instance with n = 100 qubits 
and m = 84 clauses is solved using adiabatic evolution with 
X = 14. The plot shows the entanglement entropy of a half- 
cut and the probability in solution state vs. s. 



its final value in the solution. Physically, the algorithm 
builds entanglement up to the critical point where the 
solution is singled out and, thereon, the evolution drops 
the superposition of wrong states in the register. 

This success comes at the price of a controlled loss of 
unitarity. We plot in Fig. [21 the norm in the simulation 
as a function of x i n logarithmic scale, for instances of 
14, 18, 22 and 30 qubits. The remarkable fact is that some 
observables, like the energy, appear to be very robust 
against this inaccuracy. Our simulations also allow to 
compute the decay of the x Schmidt coeficients A„ at 
any step of the computation. Close to criticality, and for 
the central cut of the system, these can be approximately 
fitted by the law log 2 (Aa ) = b + -7= + dy/a, with 
appropriate coefïicients b, c and d. 

The ultimate goal of finding the correct solution ap- 
pears also to be very robust in the simulations we have 
performed. The exact probability of success can be calcu- 
lated in 0(n poly(x)) time as well. As a symbolic exam- 
ple, our program has solved an instance with n = 100 
qubits, that is, the adiabatic evolution algorithm has 
found the correct product state out of 2 100 ~ 10 30 for 
a hard instance with m = 84 clauses and T = 2000. The 
simulation was done with a remarkable small x — 14 <C 
2 s0 = Xmax and is presented in Fig. |31 

The same robustness of evolving towards the correct 
solution is found for any number of qubits and small x- 
We have launched a search for the minimal T that solves 
samples of n-qubit hard instances in the following way: 
for a set of small vàlues of x> we try a random instance 
with an initial e.g. T = 100. If the solution is found, we 
proceed to a new instance, and if not, we restart with a 
slower adiabatic evolution e.g. T — 200. This slowing 
down of the algorithm is performed till a correct solu- 
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Figure 4: (Color online) Accumulated statistics up to n — 60 
for T m in(n) for which each instance is solved (mean and worst 
cases). Averages are performed over 200 instances for each n, 
except for n — 50,60 with respectively 199,117 instances. 
Error bars give 95 per cent of confidence level in the mean. 
The worst cases found are shown in the inset. 



tion is found and the minimal successful T m i n is stored. 
Our rcsults are shown in Fig. 0] The average over n- 
qubit instances appears to grow sub-exponentially with 
n. In fact, a quadratic fit reproduces the data for n < 22, 
consistently with the results found in 7]. The required 
times for larger n lie below the extrapolated curve. Iso- 
lat ed instances, however, may require larger times. Since 
the worst T min found depends on the interpolating path, 
finding an instance that needs a very large T m i n is no 
counterproof for the efficiency of the adiabatic algorithm, 
as alternative paths may solve the instance in a shorter 
T 0. 

In this paper we have presented simulations of quan- 
tum computation based on matrix product states that 
can be taken up to 100 qubits. The remarkable fact that 
the algorithm finds the correct solution to a large hard 
instance and the robustness in the expected energy is 
to be contrasted with the loss of unitarity inherent to 



the local truncation scheme. This drawback may well be 
ameliorated if optimal truncations are implemented. 
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